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INTRODUCTION 


of  °f  this  Pr?Sram  Package  is  to  produce  a  set 

of  relative  receiver  functions  (RRF)  for  a  given  arrav  of 

toatib^‘rv^aCh  RRF.models  a11  n«ar  receiver  contributions 
observed  seismic  waveforms;  including  crustal 

lateraTw1^01^'  ac°ustlcr  impedance  amplification  effects, 
laterally  refracted  arrivals  and  variations  in  anelastic 

azimuth  ann'-  C1farly  the  RRF’S  must  be  functions  of  back 

sourc£hrpo?nnnciaence  an?ie'  a?d  must  be  recomputed  for  each 
source  region  under  consideration. 


■  .  This  program  p-.ckage  generates  convolution  filters 
f°ri  ^Se  in  synthetic  seismogram  or  other  foreward 
modeling  applications.  The  two  major  steps  involved  in  the 

Deronin!^?11  Deconvolution  and  Minimum  Entropy 

onvolution  (MED),  and  both  steps  are  specifically 
.  to  Produce  filters  which  may  be  applied  In  a 

?efIrinr10IVSfnSe‘  In  the  trace  ^convolution  step,  a 
station  must  be  chosen  with  spectral  character 

of  thl  ttai-  nVOlU^°n  filter  (no  zeros>'  while  the  rest 
of  the  stations  m  the  array  must  appear  more  like 

U  1°n  filters  (bandlimited) .  Processing  by  MED  adds 
additional  tapers  at  the  high  frequency  end  of  the  bandpass 

o ^convolution Rf  ilters*.^  «U  °f  the  *>"’«  *»»et.riSSS' 


PROGRAM  STRUCTURE 


The  computation  of  relative  receiver  functions  is 
implemented  as  two  separate  program  packages:  DECON  and 
FMED .  DECON  performs  the  trace  deconvolution  and  outputs  an 
average  spectral  ratio  for  each  station  pair  (secondary/ 
reference).  FMED  uses  the  trace  deconvolutions 
to  estimate  relative  receiver  functions.  The  theory  behind 
the  software  is  given  by  Hart  et  al,  (1979)  and  Lundquist  et 
al,  (1980  a,b) . 

DECON  requires  a  set  of  trace  pairs,  and  the 
accumulation  and  preparation  of  the  trace  pairs  are  often 
the  most  time  consuming  portions  of  receiver  function 
estimation.  The  trace  pairs  for  different  secondary 
stations  in  the  array  need  not  be  the  same,  but  the 
reference  station  must  be  common  for  each  array.  Each  trace 
is  Fourier  transformed;  a  ratio  of  secondary  to  reference 
trace  is  taken  and  the  ratios  are  averaged  in  log  amplitude 
sense.  To  avoid  division  by  zero,  the  ratio  is  taken  only 
over  a  specified  frequency  range;  so  bandwidth  is  an 
important  consideration  in  developing  receiver  functions. 
Other  assumptions  on  the  trace  deconvolution  process  depend 
upon  the  type  of  array  used,  and  are  discussed  in  Lundquist 
et  al. ,  (1980  a). 

FMED  estimates  relative  receiver  functions  from  the 
spectral  ratios  by  finding  common  spectral  content  in  the 
set  of  spectral  ratios  for  an  array  and  ascribing  that  to 
the  reference  station.  This  is  done  by  means  of  a 
simplicity  criterion  (Wiggins,  1978)  implemented  in  the 
frequency  domain  (Lundquist  et  al,  1980  a).  The  method  is 
an  interative  technique  to  maximize  simplicity  as  a  function 
oi  a  specific  norm.  Convergence  to  a  maximum  is  a  function 
of  reference  station,  weight  functions  applied  and  the 
quality  of  the  trace  deconvolutions.  The  maximum  found  is 
not  necessarily  unique,  but  is  often  useful.  One  filter  is 
estimated  and  applied  to  the  entire  array  of  trace  decon¬ 
volutions,  so  the  receiver  functions  estimated  are  still 
relative  in  that  they  cannot  predict  absolute  amplitudes. 
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DECON  CALLLING  SEQUENCE 
DECON 

DECINS 

LINOG1 

HAGINL 

DAT  IN 

INTERP 

VNORM 

TAP 

PLOT 

COOLB 

CDATA 

DATIN 

INTERP 

VNORM 

TAP 

PLOT 

COOLB 

CDATA 

CMPR 

FPEAK 

CONJ 

COOLB 

COR1 

PLOT 

FLNFLD 

COOLB 

PLOT 


ft-  .-yrsmQ 


FMED  CALLING  SEQUENCE 
FMED 

COOLB 

GBANDP 

COOLB 

TTAPER 

COOLB 

PLOT 

COOLB 

PREMED 

COOLB 

COOLB 

CAUS 

COOLB 

COOLB 

FMDOUT 

PLOT 

PLTMED 

PLOT 


I 


t 
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Program  DECON 


This  is  the  main  program  controlling  data  input  and 
computation  and  output  of  spectral  ratio  trace  decon¬ 
volutions.  After  reading  in  preliminary  parameters, 
the  spectral  ratios  are  computed  and  averaged  in  a  major 
loop  on  event  number.  Finally,  the  result  is  output  to 
diskfile  and  plot. 

Calls:  DECINS  -  Computes  spectral  ratio  of 

secondary  station  instrument 
response  to  that  of  reference 
station. 


DAT  IN 

VNORM 

TAP 


Inputs  and  reinterpolates  one 
seismogram. 

Computes  the  transparency  norm 
of  the  input  seismogram. 

Applies  cosine  tapers  to  front 
and  back  of  seismic  trace. 


COOLB 


Utility  Fast  Fourier  Transform 
routine . 


CDATA 

CMPR 

FPEAK 

C0R1 

FLNFLD 


Finds  effective  bandwidth  of 
seismogram  spectrum. 

Computes  the  RMS  difference 
between  amplitude  spectra  of 
current  spectral  ratio  and 
previous  one. 

Filters  spectral  ratio  with 
Bartlett  filter  and  applies  a  time 
shift.  Returns  time  domain 
signal . 

Correlation  routine  to  find 
time  shift  required  to 
maximize  fit  between  two  time 
domain  sequences. 

Fills  out-of-band  spectral 
components  with  zero  and  folds 
about  Nyquist  to  complete 
spectrum. 


•  /-7;rw; 


c 
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TO,  TOS 


Reads : 


TG,  TGS 

S,  SS 
HO,  HOS 


Seismometer  free  periods  for 
reference  and  secondary  stations, 
respectively.  If  TO,  TOS  less 
than  zero,  then  TG,  TGS 
(below)  will  be  instrument 
response  number  in  a  file  of 
instrument  responses. 

Periods  of  galvanometers 
unless  TO,  TOS  less  than  zero. 

Coupling  between  TO  and  TG. 

Damping  constants  of 
seismometers . 


HG,  HGS 

- 

Damping  constants  of  galvonometers . 

XLIN ,  XLINS 

- 

Inductances . 

TDECON 

Length  of  DECON  trace  in  time 
domain  stack.  Frequency  domain 
average  of  spectral  ratio  is 
independent  of  TDECON. 

DT 

- 

Sample  interval 

NCOS 

- 

Length  of  cosine  taper  in  samples. 

SCALE 

- 

Plot  scale  factor. 

FI,  F4 

- 

Frequency  band  for  spectral  ratio 
and  frequency  domain  averaging. 

F2 ,  F3 

- 

Frequency  band  for  time  domain 
stack.  F1<F2<F3<F4 . 

ISDAC 

- 

Appropriate  to  SGI  data  input 
types,  only. 

ISP 

- 

Specifies  output  plot  requests. 

ITAPDT 

Length  of  linear  tape:  of  the 
ends  of  data  traces.  Rarely 
other  than  unity. 

ICMPS 

“ 

Flags  request  for  comparison 
between  current  spectral  ratio 
and  a  previous  run.  Rarely  used. 

START 

- 

Time  of  window  start  of  reference 
trace . 

TSTP 

- 

Stop  time  for  reference  trace  window 
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POL 

Polarity  flag.  POL  equal  to 
zero  causes  input  trace  to  be 
inverted. 

START1 

- 

Window  start  time  for 
secondary  station. 

TSTP1 

- 

Window  stop  time  for 
secondary  station. 

(NOTE:  Four  lines  of  input  are  read  in  loop  on  event 

number.  These  are  the  reference  station  window  parameters, 
the  reference  station  file  name  (read  by  DATIN),  the 
secondary  station  window  parameters  and  the  secondary 
station  file  name.  This  input  format  continues  until  START 
is  less  than  -999.) 

TSHFT  -  Time  of  circular  shift  of 

frequency  domain  averaged  traced 
deconvolution 

TZ  -  First  TZ  seconds  of  trace 

deconvolution  set  to  zero. 

TCOS  -  Next  TCOS  seconds  smoothed 

with  cosine  taper. 

I FLAG  -  Defines  convolution  option. 

<0  apply  trapezoidal  filter 
=0  apply  instrument  response 
>0  apply  Bartlett  filter. 

Though  options  are  offered, 
program  is  nearly  always  run  with 
Bartlett  filter  option. 

TWD  -  Triangle  width  in  seconds  for 

definition  of  the  Bartlett  filter. 

Outputs :  DECFFT  -  Final  frequency  domain 

averaged  trace  deconvolution 
after  time  shift  and  tapers 
applied.  This  is  input  file 
for  program  FMED . 


I 

I 


Subroutine  DECINS  (INST,  RS,  ARGS,  IF1,  IF4,  DT) 

Does  deconvolution  of  instrument  responses  to  permit 
equalization  of  instruments  during  trace  deconvolution.  The 
response  functions  are  either  computed  or  interpolated, 
depending  upon  request. 

Called  by:  DECON 


Calls : 

LINOG1 

Linear  or  logarithmic 
interpolation  routine.  Used 
when  response  function 
read  from  a  data  file. 

HAGINL 

— 

Computes  instrument  response 
using  instrument  constants  read 
by  DECON. 

Input : 

RS,  ARGS 

- 

dummy  work  space 

IF1 ,  IF4 

- 

frequency  limits  for  ratio  of 
instrument  responses 

DT 

sample  interval .  Used  with 

MAXX  (in  parameter  statement) 
to  define  DF,  the  frequency  sample 
interval . 

Output: 


INST 


the  spectral  ratio  of 
secondary  instrument  to 
reference  instrument 


Subroutine  LIN0G1(N,  X,  Y,  XZ,  Yl,  I) 


Linear 
on  input. 

Called  by: 

Calls  : 

Input  : 


Output 


interpolator,  but  assumes  X  scale  is  logarithmic 
scale  may  be  linear  or  logarithmic. 

DECINS 


none 

N 

X 

Y 

XZ 

I 


Yl 


index  of  last  point  in  X  array 

input  abscissa  array  (frequency) 

input  ordinate  array  (amplitude  or 
phase  response) 

desired  abscissa  point 

index  of  last  X  value  used  in 
previous  call  Minimizes  search  for 
X  values  bracketing  X2 . 

interpolated  ordinate  value 
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Subroutine  HAGINL  (F,  PT,  PG,  S,  HP,  HG,  XL,  R,  ARG) 

Computes  instrument  response  using  Hagiwara  formula. 


Called  by: 

DECINS 

Calls 

None 

Input  : 

F 

- 

frequency  at  which  response  desired 

PT 

- 

seismometer  period 

PG 

- 

galvanometer  period 

S 

- 

coupling  between  PT  and  PG 

HP 

- 

damping  of  galvanometer 

XL 

- 

inductance 

Output  : 

R 

- 

amplitude  response  at  frequency  F 

ARG 

- 

phase  response 

C 
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Subroutine  DATIN  (XI,  DT,  START,  NP,  IFLG) 


Inputs  one  seismogram,  reinterpolates  if  necessary  and 
subtracts  mean.  This  routine  is  highly  machine  and  user 
dependent,  and  may  be  redefined  for  each  new  format  or  data 
type  used. 

Called  by:  DECON 

Calls  :  system  routines  for  file  manipulation  and 

direct  disk  reads. 


INTERP 

linear  interpolator 

Input  : 

DT 

sample  interval 

START 

window  start  time  in  seconds 

NP 

number  of  points  in  desired 
output  array 

IFLG 

Pertinent  only  to  SGI  PRIME 
data  files.  Set  =  0 

Reads  : 

seismogram  and 
Exact  details 

descriptive  information, 
are  machine  and  user  dependent 

Output  : 

XI 

desired  seismogram 

C 
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Subroutine  INTERP  (X,  Y,  KK1 ,  XI,  DEL.  KK2 ,  XOUT) 

Linear  interpolation  subroutine. 

Called  by:  DATIN 

Calls  :  none 

Input  :  X  -  array  of  abscissa  values 

Y  -  array  of  ordinate  values 

KK1  -  length  of  X  and  Y  arrays 

in  samples 

XI  -  start  time  of  desired  output 

array  in  seconds 

DEL  -  sample  interval  of  desired 

output  array 

KK2  -  number  of  points  in  desired 

output  array 

Output  :  XOUT  -  interpolated  output  starting 

at  time  START  and  sampled  at 
time  interval  DT 


Subroutine  VNORM  (X,  V,  NP) 


Computes 

Called  by: 
Calls  : 
Input  : 

Output  : 


the  transparency  norm  defined  by: 

4 

V  = 

(1YT) 

DECON 

none 

X  -  time  domain  array 

NP  -  length  of  X  in  samples 

V  -  transparency  norm 
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Subroutine  TAP  (X,  ISTP,  NCOS,  ITAPDT ) 


Applies  cosine  tapers  to  front  and  back  of  a  time 
domain  array  and  linear  tapers  to  the  back  of  the  array. 


Called  by: 

DECON 

Calls  : 

None 

Input  : 

X 

- 

array  to  be  tapered 

ISTP 

- 

length  of  X  in  samples 

NCOS 

- 

length  of  cosine  tapers 

ITAPDT 

- 

length  of  linear  taper 

Output: 

X 

- 

tapered  version 

14 


Subroutine  (PLOT) 

System  dependent  time  domain  plot  routine. 


Called 

by:  Various 

Calls 

:  System  plot 

routines 

Input 

:  X 

array  to  be  plotted 

DT 

sample  interval 

N 

length  of  X  in  samples 

SCALE 

scale  value  such  as  number 
of  mm/sec. 

Output 

:  hard  copy  plot 

15 
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Subroutine  COOLB  (NN,  XX,  SIGNI) 


Utility 
Transforms . 

Called  by: 

Calls  : 

Input 


Output 


subroutine  for  computation  of  Fast  Fourier 


Various 

None 

XX 

NN 

SIGNI 

XX 


complex  array  to  be 
transformed 

power  of  2  describing  length  of  XX 

Forward  transform  done  when 
SIGN!  =  -1.0; 

inw  transform  done  when 
Oft-  •;  i.0. 

•:r  s *>r  ;'-::*med  version 
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Subroutine  CDATA  (C,  I FI,  IF4,  WTL,  XMAX) 


Searches  amplitude  spectrum  for  maximum  value  of 
spectrum  and  for  frequencies  at  which  spectrum  has  decayed 
to  the  water  level  relative  to  the  maximum. 


Called  by: 

DECON 

Calls 

None 

Input  : 

C 

I  FI 

IF4 

WTL 

Output  : 

I  FI 

IF4 


XMAX 


complex  spectrum  to  be 
evaluated 

index  of  minimum  frequency  to  be 
tested 

index  of  maximum  frequency  to  be 
tested 

water  level 

either  as  input  or  index  of 
frequency  at  which  water  level 
is  reached  on  low  frequency 
side  of  maximum  value. 

either  as  input  or  index  of 
frequency  at  which  water  level 
is  reached  on  high  frequency  side 
of  maximum  value. 

maximum  spectral  amplitude 
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Subroutine  CMPR  (I FI,  IF4) 


Compares  amplitudes  of  two  spectra  input  in  common  and 
outputs  RMS  difference  between  the  two  to  a  printer  file. 

Called  by:  DECON 

Calls  :  None 

Input  :  1F1,  IF4  -  frequency  limits  for 

comparison 


Output 


None 


Subroutine  FPEAK  (PEAK,  I FI,  IF4,  DT,  SHFT) 

Applies  a  Bartlett  filter  to  a  spectrum  input  in 
common,  applies  a  time  shift,  inverse  transforms  and  returns 
time  domain  array. 


Called  by: 

DECON 

Calls  : 

CONJ 

- 

folds  complex  spectrum  about 
the  Nyquist  frequency 

COOLB 

- 

utility  Fast  Fourier  Transform 
routine 

Input  :  I FI , 

IF4 

- 

frequency  limits 

DT 

- 

sample  interval 

SHFT 

- 

desired  time  shift  in  seconds 

Output  : 

PEAK 

- 

filtered  array 

19 


Subroutine  CONJ  (C,  LX1 ) 

Folds  complex  spectrum  about  Nyquist  frequency  to 
obtain  complete  frequency  domain  representation. 

Called  by:  FPEAK 

Calls  :  None 

Input  :  C  -  complex  spectrum  defined  from 

D.C.  to  the  Nyquist  frequency 

LX1  -  index  of  the  Nyquist  frequency 

Output  :  C  -  complete  spectrum 
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Subroutine  C0R1  (PEAK,  PEAK1,  JSHFT,  NPTS,  IMAX) 

Routine  to  estimate  the  shift  required  to  maximize  fit 
between  two  time  domain  arrays. 


Called  by: 

DECON 

Calls  : 

None 

Input  : 

PEAK, 
PEAK1 , 

- 

the  time  domain  arrays  to  be 
compared 

JSHFT 

- 

maximum  shift 

NPTS 

- 

number  of  points  to  compare 

Output  : 

IMAX 

- 

shift  value  at  maximum  correlation 

Subroutine  FLNFLD  (C,  IF1,  IF4) 


Routine  to  fill  out-of-band  components  of  a  bandlimited 
spectrum  with  zeros,  then  fold  about  Nyquist  frequency. 


Called  by: 

Various 

Calls  : 

None 

Input  : 

C 

complex  spectrum  to  t 
completed 

IF1 ,  IF4  - 

indices  of  bandlimit 
frequencies 

Output  : 

C 

complete  spectrum 

Program  FMED 

This  is  the  main  program  controlling  input  of  trace 
deconvolutions  (from  DECON)  and  computation  of  relative 
receiver  functions.  The  program  reads  in  DECFFT  files  for 
each  secondary  station,  initializes  a  corresponding  array 
for  the  reference  station,  and  then  attempts  to  maximize  the 
simplicity,  or  minimize  the  entropy,  of  the  set  of  trace 
deconvolutions  by  maximizing  the  varimax  norm.  Effectively 
this  determines  common  spectral  character  in  the  array  of 
DECON  transfer  functions,  and  puts  that  spectral  content 
into  the  reference  station  receiver  function.  Since  one 
filter  is  determined  and  applied  to  all  input  traces,  the 
resulting  receiver  functions  are  still  relative  in  that  any 
spectral  content  common  to  all  seismograms  will  still  be 
missing.  In  particular,  absolute  amplitude  is  not  predicted 
by  these  filters. 


Calls  :  COOLB 

GBANDP 

TTAPER 

PLOT 

PREMED 

I 

Reads  :  NSTA 

LSG 

NITER 

DT 

( 


utility  Fast  Fourier  Transform 
routine 

Gaussian  bandpass  routine 

applies  several  time  domain 
tapers 

system  dependent  plot  routine 

performs  the  iterative 
computation  and  application 
of  the  minimum  entropy  decon¬ 
volution  filter. 

number  of  stations  in  the 
array.  Includes 
reference  station. 

number  of  points  output  for  ' 
each  station  at  each  iteration 

number  of  iterations 

sample  interval 

initialize  filter  array  as  a 
spike  at  sample  IFP 


IFP 


* 


I 


I  TP 


I  REF 


TCOS 

TZ 

TEXP 


IWTOPT 

WTL 


FI,  F2 , 
F3 ,  F4, 


ISP 


ISOUT 

W 


initialize  reference  station 
trace  as  spike  at  ITP.  Reference 
trace  will  be  filtered  by  the 
same  bandpass  used  on  other 
traces.  That  is,  the  reference 
trace  is  initialized  as  the 
bandpass  impulse  response  centered 
at  sample  ITP. 

index  of  the  station  to  be  used 
as  the  reference  for  this  FMED 
run.  If  same  as  DECON  reference, 
then  IREF=NSTA.  Otherwise,  all 
traces  multiplied  by  the  inverse  of 
trace  IREF,  effectively  changing 
reference  station. 

time  of  beginning  of  cosine 
taper  of  end  of  traces 

cosine  taper  ends  and  zero  trace 
begins  at  time  TZ 

generate  causal  time  domain 
signal  by  zeroing  trace  up  to 
arrival  time.  TTAPER  searces 
backward  from  TEXP  to  a  zero 
crossing  and  smoothly  zeros 
all  earlier  points. 

option  number  for  definition  of 
weight  function 

percent  maximum  spectral 
amplitude  for  prewhitening 


frequencies  for  bandpass  filter 

specifies  stations  for  which 
spectra  will  be  plotted 

GT.O  plot  both  input  and  final 
output 

EQ.O  no  plots 

LT.O  plot  final  output  only 

NE.O  time  domain  receiver 
functions  output  to  diskfiles 

weight  functions.  Must  read, 
even  though  a  different  weight 
option  is  used. 
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Input 

File 

names  :  Usually  DECFFT  files  from  DECON. 

SCALE  -  Gain  factor  (Z  =  Z/scale).  May 

include,  if  required,  correction 
for  geometric  attenuation. 

Output 

file 

names  :  Used  for  plot  labels  even  if 

ISOUT  =  0. 

Outputs:  Plots  of  spectra,  if  requested. 

All  other  output  is  done 
by  PREMED. 
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Subroutine  GBANDP  (X,  KF1 ,  KF2 ,  KF3 ,  KF4 ) 


Applies  a  simple  Gaussian  bandpass  filter  to  the  input 
complex  array. 

Called  by:  FMED 

Calls  :  None 

Input  :  X  -  array  to  be  filtered.  Length 

defined  by  parameter  statement 
to  be  consistent  with  FMED. 

KF1 ,  KF2  -  indices  of  frequencies  with 
zero  amplitude  and  unit 
relative  amplitude, 
respectively,  or  the  low 
frequency  end  of  the  bandpass 

KF3 ,  KF4  -  indices  of  frequencies  with 
unit  relative  amplitudes  and 
zero  amplitude,  respectively, 
on  the  high  frequency  end  of 
the  bandpass 

Outputs  :  x  -  filtered  version 


Subroutine  TTAPER  (Z,  TCOS,  TZ,  TEXP,  DT) 


Routine  applies  several  time  domain  tapers  to  a  complex 
time  domain  array.  (1)  Zeros  all  imaginary  parts.  (2) 
Applies  cosine  taper  from  time  TCOS  to  TZ  at  end  of  trace. 
(3)  Zeros  all  amplitude  for  times  greater  than  TZ.  (4) 
Applies  an  exponential  taper  to  front  of  trace  from  zero 
crossing  before  TEXP  to  first  sample. 

Called  by:  FMED 

Calls  :  None 


Input  :  Z 


complex  time  domain  array  to 
be  filtered 


TCOS 

TZ 

TEXP 


DT 


time  to  start  cosine  taper 

time  to  end  cosine  taper  and 
begin  zero  fill 

Routine  searches  back  in  time 
from  TEXP  to  find  a  zero 
crossing.  Trace  amplitudes 
before  the  zero  crossing  are 
smoothed  exponentially. 
Primary  use  is  for  causal 
truncations.  Rarely  used  in 
practice . 

sample  interval 


Output  :  Z 


tapered  version 
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Subroutine  PREMED  (IREF) 


Routine  to  compute  filter  which  maximizes  the  verimax 
norm  in  the  frequency  domain.  The  theory  and  a  discussion 
of  this  application  are  given  by  Wiggins  (1978)  and 
Lundquist  et  al . ,  (1980  a,b).  The  filter  is  initialized  in 
FMED  as  a  delayed  delta  function.  The  filter  is  applied  in 
this  routine  to  each  DECON  spectral  ratio  and  the  varimax 
norm  is  computed.  The  filter  is  then  reestimated  based  upon 
the  norm  and  reapplied  to  the  data.  The  process  is  done  for 
a  specified  number  of  iterations. 

Note  that  the  convergence  of  the  technique  is  not 
guaranteed.  It  is  not  sufficient  to  declare  a  minimum 
change  in  the  norm  and  let  the  program  iterate  until  that 
threshold  is  found.  Personal  attention  is  required  to 
maintain  adequate  control  of  the  process. 


Called  by: 

FMED 

Calls  ; 

COOLB 

- 

utility  Fast  Fourier  Transform 
routine 

CAUS 

- 

applies  causal  tapers  to  the 
filter 

FMDOUT 

outputs  relative  receiver 
functions  and  descriptive 
information  to  diskfile 

Input  : 

IREF 

- 

trace  number  used  for 
reference  trace 

All  flags  and  parameters  read  by  FMED  are  made 
available  to  PREMED  in  common  blocks. 

Output  :  Relative  Receiver  Functions  are  output  to 

diskfiles,  if  requested,  at  the  last  iteration. 

A  disk  file  is  generated  with  the  output  from 
each  it  iteration.  This  will  be  used  for 
input  to  PLTMED  for  plotter  output. 

A  separate  diskfile  is  generated  which  allows 
monitoring  of  the  varimax  norms. 
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Subroutine  CAUS  (F,  NN) 


Routine  tapers  the  FMED  filter.  Filter  is  transformed 
to  time  domain;  all  samples  before  first  zero  crossing  are 
set  to  zero;  all  samples  in  the  second  half  of  the  array  are 
set  to  zero,  and  the  array  is  transformed  back  to  frequency. 

Called  by: 

PREMED 

Calls 

COOLB 

- 

utility  Fast  Fourier  Transform 
routine 

Input  : 

F 

- 

FMED  filter  in  complex 
frequency  domain 

NN 

- 

power  of  2  defining  length  of 

F 

Output  : 

F 

_ 

tapered  version 

Subroutine  FMDOUT  (II,  NOUT)) 

Routine  to  output  a  relative  receiver  function  and  its 
descriptive  information  to  a  diskfile. 

Called  by:  PREMED 

Calls  :  System  dependent  file  open  and  close  routines. 

Input  II  -  station  number  for  current  RRF 

output 

NOUT  -  number  of  points  to  output 

All  flags  and  parameters  read  by  FMED  are  made 
available  to  FMDOUT  in  common  blocks, 
including  name  of  output  files. 

Output  :  Diskfiles  containing  NOUT  points  of  relative 

receiver  function  for  station  II,  along 
with  descriptive  information. 
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Program  PLTMED 

This  program  plots  the  results  of  FMED.  The  "  .suits 
for  each  station  at  each  iteration  were  written  to  .ile 
MEDPR  by  PREMED,  so  the  user  may  plot  some  or  all  of  the 
output.  In  normal  usage,  the  analyst  would  first  examine 
file  FMEDFR  which  contains  the  varimax  norms  for  each 
station  and  each  iteration.  Convergence  and  instablility  are 
indicated  by  those  norms,  and  the  amount  of  plotting  can  be 
tailored  to  the  amount  of  change  in  the  waveform  as 
determined  by  change  in  the  norms. 

Calls  :  System  plot  routines 

Reads  :  IPITER  -  Flag  determining  which 

iterations  to  plot. 

IPSTA  -  Flag  determining  which 

stations  to  plot. 

Output  will  be  generated  only 
when  both  IPITER  and  IPSTA  are 
nonzero. 

MEDPR  _  disk  file  containing  output  of 

FMED  for  each  station  at  each 
iteration.  The  input  data  is 
present  as  iteration  zero. 

Output  :  Hard  copy  plots. 
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